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ABSTRAGT 

We carry out direct numerical simulations of turbulent astrophysical media that explicitly track 
ionizations, recombinations, and species-by-species radiative cooling. The simulations assume solar 
composition and follows the evolution of hydrogen, helium, carbon, oxygen, sodium, and magnesium, 
but they do not include the presence of an ionizing background. In this case, the medium reaches 
a global steady state that is purely a function of the one-dimensional turbulent velocity dispersion, 

(JiD, and the product of the mean density and the driving scale of turbulence, nL. Our simulations 
span a grid of models with ctid ranging from 6 to 58 km s“^ and nL ranging from 10^^ to 10^^ cm“^, 
which correspond to turbulent Mach numbers from M = 0.2 to 10.6. The species abundances are 
well described by single-temperature estimates whenever M is small, but local equilibrium models 
can not accurately predict the global equilibrium abundances when M > 1. To allow future studies 
to account for nonequilibrium effects in turbulent media, we gather our results into a series of tables, 
which we will extend in the future to encompass a wider range of elements, compositions, and ionizing 
processes. 

Subject headings: ISM: abundances, ISM: atoms, astrochemistry, turbulence 


1. INTRODUCTION 


Turbulence is ubiquitous in astrophysics, where the 
Reynolds number (Re), the ratio of the inertial forces to 
the viscous forces, is often orders of magnitudes higher 
than found on the Earth. In the intergalactic medium, 
for example. Re typically exceeds 10^, and in the warm 


ionized inters tellar medium Re > 10^ ( [Bragins^ 


1958 


Spitzer 1962), whereas the transition from l ammar to 
turbulent ff ow occurs at Re ~ 3 x 10^ (e.g., Orszag & 
Kells||198Q ). In addition, in many astrophysical regimes. 


rapid cooling causes turbulent velocities to exceed the 
sound speed, and within such supersonic turbulence, ran¬ 
dom motions can compr ess a fraction of the material to 
very high densities (e.g., Padoan et al.|1997 Mac Low & 


Klessen 
plex, miilti-p. 


20041 Eederrath et al. 2008|), producing a com- 
Llti-pnase medium. 


Eurther complicating the picture is the fact that the 
recombination and collisional ionization times for many 
species are long with respect to the “eddy turn-over 
time” on which existing turbulent motions decay and new 
turbulent motions are added. Thus, the conditions expe¬ 
rienced by a parcel of gas may change before equilibrium 
can be reached, such that the ionization structure of the 
medium will depend not only on the small scale density 
and temperature distribution, but on the velocity distri¬ 
bution as well. Eor these reasons, the turbulent structure 
of gas can have a significant impact on line emission and 
absorption diagnostics. Most interpretations of observed 
spectra, however, do not take into account these multi¬ 
phase and non-equilibrium effects, and full galaxy simu¬ 
lations can not resolve the relevant small-scale structures 
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to include them. 

Eurther more, there are several reasons to expect that 
nonequilibrium effects are important in interpreting cur¬ 
rent observations. Eor example, emission line studies of 
high star-formation rate, ultraluminous infrared galaxies 
(ULIRGs) have progressed to the point that a number 
of sens itive line diagnostics can now be considered. Re¬ 
cently, I Soto & Martin (|2012 ) were able to measure emis¬ 
sion lines from material in and around 39 ULIRGs, using 
the measurements to construct line ratios that are rela¬ 
tively insensitive to the presence of dust, but hig hly con¬ 


straining of the sources of ph otoionization (e.g., Kewley 


et al.||2006 

Allen et al. 

2008). They found that material 

up to ~ lU kpc from 

bhe galaxy centers was primarily 


heated by shocks rather than photoionization, as would 
be expected for strong ly turbulent media . 


At larger distances, Werk et al. 


used the Cos¬ 


mic Origins Spectrograph (jGreen et aL|| 26 l 2 ) to measure 
low and intermediate ionization state ions m the circum- 
galactic medium (G GM) within ^ 100 kpc o f low-redshift 
galaxies (see also, 


Tumlinson et al 


2013 Peeples et al.|2014[ ). To model t 


1[2M3| 

heseabi 


Werk et al 

___ sorbers they 

adopted several important assumptions, namely: (i) that 
the ions were co-spatial and arose from a single phase; 

(ii) that the medium was in ionization equilibrium; and 

(iii) that the medium was primarily photoionized. Sur¬ 
prisingly, to match the observations with such assump¬ 
tions, the models had to adopt large ionization param¬ 
eters, defined as the ratio of ionizing photon density to 
the hydrogen density. Given the observed range of meta- 
galactic and host galaxy fiuxes, these ratios corresponded 
to densities and pressures over two orders of magnitude 
lower than expected from hydrostatic balance. However, 
given the long recombination and cooling times in the 
diffuse GGM, even moderate energy input from decaying 
turbulence may be sufficient to substantially change this 
picture. 

It is with these issues in mind that we have carried out 
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a comprehensive numerical survey of the atomic struc¬ 
ture and observational properties of turbulent astrophys- 
ical media. While a large number of high-resolution 
studies of turbulence exist in the literature, the major¬ 
ity of these are incompressible simulations ca rried out in 
the context of fluid-dynamics research ( e .g., [Vincent fc 
Meneguzz^ 1991| |Moin fc Mahesh 1998[ Ishihara et at 
2QQ9p , and compressible, isothermal siuiulations carried 


out m the context of interstellar medium research (e.g.. 


Kritsuk et al.||2007l IFederrath et al.||2008l 

Schmidt et al. 

2009 

Pan fc Scannapieco||2010t Sur et al. 

2014||. Walch 

et al.]i 

201 ip carried out simulations of both continuously- 


0.001 solar metallicity material in a medium with a one¬ 
dimensional turbulent velocity dispersion, ctid ~ 30 in a 
500 parsec cubical box with a mean h ydrogen density of 
1 cm“^, using the chemical network of Glover fc Jappsen 


(2007) to track the nonequilibrium c hemistry associated 
with Hydrogen an d helium (see also Walch et al.||2014 ). 


Saury et al. (2014), on the other hand, studied the struc- 
ture of therniially bistable continuously-driven turbulent 
gas, using a cooling rate given as n^A(T) and a heating 
rate given as nr(T), where n was the local number den¬ 
sity and A and F were global functions accounting for 
atomic cooling, recombination on interstell ar grains, and 
photo -electric heating of small dust grains (Wolfire et al. 
200 ^. 

Here we carry out exact calculations of atomically- 
cooled astrophysical media that explicitly track 
continuously-driven turbulent motions, radiative cooling 
by atomic species, and the nonequilibrium ionization 
and recombination of several elements, for a grid of 
solar-metallicity models with ctid ranging from 6 to 
58 km s“^ and the product of the mean density and 
the turbulent driving scale ranging from 10^^ to 10^^ 
cm“^. In this first paper in this series, we track the 
nonequilbrium ionization structure of hydrogen, helium, 
carbon, oxygen, sodium, and magnesium, focusing on 
the limit in which ionization is purely collisional. Taking 
advantage of the scaling properties of this case, we 
are able to fully span the relevant range of conditions 
experienced in atomically cooled astrophysical plasmas: 
cataloging their physical properties for comparison with 
more idealized simulations and tabulating their species 
mass fractions and Doppler parameters for use in the 
interpretation of future theoretical and observational 
studies. 

The structure of the paper is as follows. In §2 we de¬ 
scribe our atomic chemistry and cooling routines and 
their associated tests. In §3 we present our simulation 
setup and initial conditions, and in §4 we describe our 
results, taking particular note of their probability den¬ 
sity functions and the effect of resolution and thermal 
conduction. Concluding remarks are given in §5. 

2. NUMERICAL METHOD 

All simulations were p erformed with FLASH version 
4.0.1 ( Fryxell et al.||2QQQ| ), a publicly-available hydrody- 
namics code. T'o ensure tne stability of the code as turbu¬ 
lence develops, we employ a hybrid Riemann solver which 
uses both an extremely accurate but somewhat f ragile 
Harten-Lax - van Leer-Co ntact (HLLC) solver (e.g., Toro 


feldt et al. 1991). The HLLC solver is a modification 
to the HLL solver that includes the missing shear and 
contact waves, and it produces solutions that most ac¬ 
curately capture contact discontinuities. However, in 
regions with strong shocks or rarefactions, the HLLC 
solver can fail, and in such situations, we switch to 
the positivity-preserving HLL solver. Magnetohydro dy¬ 
namic effects were not included in this study. 

In order to accurately determine the atomic proper¬ 
ties of turbulent media, we also added two new capabil¬ 
ities to the code: a non-equilibrium ionization package 
that tracks the ionization state of several atomic species, 
and a cooling routine that takes into account the cooling 
from each of these ionized states individually. In this sec¬ 
tion we describe our numerical implementation of each 
of these new capabilities. 

2.1. Atomic Ionization 

Our ionization network tracks the impact of 36 sepa¬ 
rate reactions acting on 24 species and 6 elements: hydro¬ 
gen (H, H+), helium (He, He+, He^+), carbon (C-C^+), 
oxygen (0-0^+), sodium (Na, Na+), magnesium (Mg, 
Mg+, Mg^+), and free electrons (e“). For each species, 
we track several reactions, including radiative recombi¬ 
nations, dielectronic recombinations, and ionization due 
to electron impacts. A summary of the reactions consid¬ 
ered as well as their source is given in Table 

Our implementation of this network follows the over¬ 


all approach we have adopted in previous studi es ( Gra; 
fc: Scannapiecol ^1Q| [Scannapieco et ar]|2Q12| [Cray 
Scannapieco||2013|). Each species is labeled by an index 
A such that species i has protons, A^ nucleons, and 
a mass density of pi. We then define a mass fraction of 
species i as = Pi/p, where p = ^ipi^ and define the 
molar mass fraction as = X^/A^. For each species, a 
continuity equation for the molar mass fraction is then 
given as Yi = Ri, where Ri is the total reaction rate due 
to all reactions. 

Because of the complex ways that reaction rates de¬ 
pend on temperature and because of the large range of 
possible abundances, the resulting rate equations are of¬ 
ten ‘stiff,’ z.e., the change in timescales can be very dif¬ 
ferent from one species to another. This requires im¬ 
plicit or semi-implicit methods to track all the relevant 
reactions throughout our simulations. To handle this, 
we use a variable-order B ader-Deufihard method (e.g., 
Bader fc Deuflhard||1983|) whic h is well suited to prob¬ 
lems with d imension > lU (e.g., I [Press et al.|1992 ||Bovino 


et al 


et al. 1994 Toro 1999) and a more robust, but more 


dittusive Harten Lax and van Leer (HLL) solver (Ein- 


2013). The network presented here also has the 
advantage of being sparse. That is, the Jacobian, de¬ 
fined as Jij = dYijdYj, is mostly comprised of zeros, 
with nonzero values falling on or near the diagonal. This 
allows us to use the M A28 sparse linear algebra package 
included with FLASH (iDuff et al.||1986|) to compute the 
matrix inverses required by the Bader-Deufihard scheme, 
leading to a very fast and efficient implementation. In 
general, the chemistry package runs slightly faster than 
the hydrodynamics. 

As the species evolve over a given time step, the tem¬ 
perature can change as the internal energy changes from 
ionizations and recombinations. Since the reaction rates 
are strong functions of temperature, the network can be¬ 
come unstable if too large a time step is used, especially 
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in regions of strong cooling. To ensure testability but 
allow the simulation to evolve at the hydrodynamic time 
step, we sub cycle the rate equations within each cell on 
a network time step, defined as: 


tnet = min a 


Yi+0.1Y+ \ 




(1) 


where a is a constant that controls the maximum abun¬ 
dance change allowed, which we default to 0.1. Using the 
current species abundances, the instantaneous change in 
abundances U is computed using the analytic ordinary 
differential equations at the current temperature. Note 
that a small fraction of is added as a buffer to pre¬ 
vent species with very small, but rapidly changing, abun¬ 
dances from causing prohibitively small time steps and 
excessive subcycling. In addition to this network sub- 
cycling, the Bader-Deufihard method includes its own 
internal sub cycling and time step controls. 

At the beginning of each cycle, tnet is compared to 
the hydrodynamic time step. If the hydrodynamic time 
step is smaller, then no subcycling is done and the hy¬ 
drodynamic time step is used to update the molar mass 
fractions. On the other hand, if the network time step is 
smaller, then the molar mass fractions are updated using 
the tnet: which is then recalculated and compared to the 
remaining hydrodynamic time step. This cycle continues 
until the network has advanced for the full hydrodynamic 
time step. 

To test the implementation of our network, we emu¬ 
lated a series of equilibrium models by fixing the density 
and temperature in the simulation and running each case 
until the species reached equilibrium. We ran a series of 
seven such models that spanned the temperature range 
between 10^ and 10^ K, at a fixed density of 2.0x10“^^ 
g cm“^. All species were assumed to be neutral at the 
start of each simulation, and the final abundances were 
compar ed to equilibrium m odels from Cloudy (version 
lO.OI) (jFerland et al.||I998|), using the “coronal equilib¬ 
rium” command that enforces only collisional ionizations. 

We found that we matched the results from Cloudy 
very well over a wide range in temperature for all species. 
However, we did not match precisely for certain higher 
ionization states, as seen, for example, in the middle 
panel of Fig. but this is to be expected since our 
network does not follow every ionization state for some 
elements. We did find excellent agreement for those ele¬ 
ments for which the ionization state changes by orders of 
magnitude with a modest increase in temperature. For 
the purposes of comparison, we summed together the 
higher ionization states that we do not follow and as¬ 
signed their abundance to the highest state that we do 
follow. For example, we track the abundance of C^+ and 
not the two higher states. Therefore, the (dashed) navy 
line in the middle left panel of Fig. [T] represents the sum¬ 
mation of C^+, C^+, and C^+. 

2.2. Cooling 

The second capability we implemented is individual 
cooling rates associated with each ion in the chemical 
ne twork. To do this , we fol lowed the procedure presented 
in Gnat & Ferland| (2012) who compiled the ion-by-ion 
cooling efiiciencies for several atomic species for between 
10^-10^ K. These efficiencies include all the cooling pro- 


cesses cons i dered in Cloudy, as described in Osterbrock 


& Ferland (2006). This includes collisional excitations 
followed by line emission, recombinations with ions, col¬ 
lisional ionizations, and thermal bremsstrahlung. Here 
we briefly outline this procedure, which we used to ex¬ 
tend their work down to 5x10^ K for the ions in our 
network. Again, making use of Cloudy, we construct a 
grid of models that span 5x10^ to 10^ K in tempera¬ 
ture for each ion i of element E. Each model is then 
composed of only hydrogen, electrons, and the ion under 
consideration. To ensure cooling from hydrogen is sup¬ 
pressed, the hydrogen density is set to 10^^ times smalle r 
than the ion density. Following Gnat & Ferland (2012), 
the elemental ion and electron density is set to = 
rig- = I.O cm“^ and uh = 10“^^ cm“^. This produced 
a cooling rate for each ion as a function of temperature, 
and a comp arison of these rates to those of |Gnat fc: Fer-| 


land (2012) show identical results where the temperature 
range overlapped. 

For those species in which we do not follow every ion¬ 
ization state in our simulation individually (such as C^+, 
which captures the joint abundance of C^+, C^+, and 
C^+) we formed a composite cooling curve as 


A(T) 


E,n,(r)V(r) 

EMt) 


( 2 ) 


where Ei(T) is the cooling rate of ion i and ni{T) is the 
relative abundance of ion i in equilibrium at a tempera¬ 
ture T. This composite is then used as the cooling rate 
for the highest ionization state followed. Similarly, we 
included cooling from nitrogen, neon, silicon, sulphur, 
calcium, and iron. Again, we assumed that the relative 
ionization abundances were determined by the collisional 
ionization equilibrium. This final cooling curve has the 
form 


Aother(F) = 

3 


E,n,-,(T)A,-,(T) 


(3) 


where j loops over nitrogen, neon, silicon, sulphur, cal¬ 
cium, and iron. 

With the addition of cooling, an important timescale is 
created, namely one that relates the total internal energy 
to energy loss rate per time as 

icooi = (4) 

Ei 


where Ei is the internal energy, Ei is the energy loss rate, 
and ac is a constant between 0 and I. As is the case 
with reaction rates, cooling rates are strong functions 
of temperature and species abundances, both of which 
can change drastically over a single chemistry time step. 
Therefore, we employ a method of subcycling the cooling 
on a cooling timescale within the chemistry routine. At 
the beginning of each cooling cycle, we calculate the cool¬ 
ing timescale using a given o^c(=0.1). If this is shorter 
than the chemistry time step, we cycle over the cooling 
timescale, recalculating the temperature, cooling rate, 
and timescale at the end of each cycle. This proceeds 
until we have reached the chemistry time step. 

Since we only consider two body interactions, the cool¬ 
ing time should scale linearly with density. To test this 
scaling, along with the overall dependence of our rou- 
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Fig. 1. — Comparison of the species abundances between Cloudy and FLASH. Each panel shows the results for a different set of atomic 
species. The lines corresponds to the Cloudy results, while the points are from FLASH. We plot the equilibrium temperature along the 
x-axis and the fractional abundance of each species, i.e. Fi = nijns where Us is the elemental abundance, for each species as the y-axis. A 
universal legend is given at the top of the figure and the same ionization state is given by the same color and line style in each panel. Top 
Left: Hydrogen. Top Right: Helium. Middle Left: Carbon. Middle Right: Oxygen. Bottom Left: Sodium. Bottom Right: Magnesium. 


tines on temperature, we ran several tests in which the 
medium was initially completely ionized and the initial 
temperature was set to 1 x 10^ K, each of which had a dif¬ 
ferent density between 5.Ox 10“^^ and 1.0x10“^^ gcm“^. 
These tests recovered the expected scaling with density, 
and also confirmed that cooling was always efficient at 
high temperatures and much less efficient at ^10^ K, 
when most elements become neutral. Finally, we per¬ 
formed additional tests to ensure that the interpolation 
was smooth and always reproduced the table values un¬ 
der the correct conditions. 

3. MODEL FRAMEWORK & INITIAL CONDITIONS 

With these procedures in place, we carried out a suite 
of simulations of turbulent media under a wide variety of 
conditions. To drive turbulence, we made use of an ar¬ 
tificial forcing term F, incorporated into the momentum 
equation as 

^ + V(pvv) + VP = pF, (5) 

where p is the density, P is the pressure, and v is the 
velocity. The forcing term was modeled as a stochas¬ 


tic Ornstein-Uhlenbeck pr ocess dUhlenbeck fc Ornstein 


19301 ISchmidt et al.||2QQ9| |Federrath et aL||2010[ [Fanlz 


Scannapieco 2Q1Q| ) with a user-dehned forcing correla¬ 


tion time tf. For all the simulations presented here, tur¬ 
bulence was driven solely though solenoidal modes {i.e. 
V-F = 0) in the range of wavenumbers 1 < I/box|k|/27r < 
3, such that the average forcing wavenumber was ^ — 
2Tbox/27r, with Lbox the size of our turbulent box, which 
was fixed at 100 parsecs on a side. This turbulence was 
always continuously driven throughout the simulation 
runtime. 

All our simulations were performed using the multi¬ 
species extension for the ideal gas equation of state, 
which calculates the important thermodynamic quan- 
tities based on the p roperties of the included species 
(Colella & Glaz||1985). In particular, the average atomic 
mass can change dramatically as the gas either recom¬ 
bines or becomes ionized. While FLASH has an adap¬ 
tive mesh capability, this is capability is not useful in a 
case such as ours in which structures of interest are dis¬ 
tributed throughout the full simulation volume. Thus all 
simulations were performed using a cubic uniform grid 
with periodic boundary conditions. 
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TABLE 1 

List of the reactions considered and their source. 


Number 

Reaction 

Source 

0001 

He++ e“^ He 

1 

0002 

C++ e-^ C 

5 

0003 

C^+T e“^ C+ 

4 

0004 

C3++ e-^ C2+ 

3 

0005 

C4++ e-^. C3+ 

2 

0006 

0++ e-^ O 

7 

0007 

O^+T e“^ 0+ 

6 

0008 

03++ e-^ 02+ 

5 

0009 

0^++ e“^ 0^+ 

4 

0010 

0^++ e“^ 0^+ 

3 

0011 

0®++ e“^ 0^+ 

2 

0012 

0'^++ e“^ 0®+ 

1 

0013 

Mg++ e“—)■ Mg 

8 

0014 

Mg2++ e“^ Mg+ 

9 

0015 

Na++ e“—)■ Na 

8 

0016 

H++ e-^ H 

12 

0017 

H + e-^ H++ 2 e- 

11 

0018 

He + e“—He++ 2 e“ 

11 

0019 

He++ e-^ He2++ 2 e" 

11 

0020 

C + e —y C++ 2 e 

11 

0021 

C++ e-^ C2++ 2 e- 

11 

0022 

C2++ e“^ C^++ 2 e“ 

11 

0023 

C3++ e-^ C4++ 2 e- 

11 

0024 

O + e-^ 0++ 2 e- 

11 

0025 

0++ e“^ 02++ 2 e“ 

11 

0026 

02++ e“^ 0^++ 2 e“ 

11 

0027 

0^++ e“^ 0^++ 2 e“ 

11 

0028 

0^++ e-^ 0^++ 2 e- 

11 

0029 

0^++ e“^ 0®++ 2 e“ 

11 

0030 

0®++ e“^ 0'^++ 2 e“ 

11 

0031 

Mg + e“—Mg++ 2 e“ 

11 

0032 

Mg++ e“^ Mg2++ 2 e“ 

11 

0033 

Na + e“—Na++ 2 e“ 

11 

0034 

He2++ e-^ He+ 

12 

0035 

He++ H ^ He + H+ 

12 

0036 

He + H+^ He++ H 

12 


nearly reach collisional ionization equilibrium. Each 
model was then run normally until it reached a global 
steady state in terms of both the hydrodynamic variables 
and the chemical abundances. 

In reactions that involve free electrons recombining 
with ions, the optical depth of the environment becomes 
im portant. If the e nvironment is optically thin (Case 
A; )Osterbrock|jl989|), the ionizing photons were allowed 
to escape, while in the optically thick case (Case B), 
the photons were reabsorbed by a nearby neutral atom, 
which has the effect of lowering overall recombination 
rate. We included these effects for hydrogen, which has 
the highest number density and provides most of the free 
electrons. To best estimate the appropriate case for each 
run, at each time step we calculated the optical depth 
as, 

Th = XHpCTuLho^/mH, (6) 

where Xh is the global neutral hydrogen mass fraction, p 
is the mean ambient density, = 3.3 x 10“^^ cm^ is the 
photoionization cross section for hydrogen, and mn is the 
mass of hydrogen. The recombination rate for hydrogen 
is then 

Kec = e-^^kA + (1.0 - (7) 


where k^ec is the new recombination rate and k^ and k^ 
are the Case A and B recombination rates respectively, 
which differ by a factor < 2. When the optical depth is 
low, 1.0 and we defaulted to the Case A rate. 

Conversely when the optical depth is large, e“^^ ~ 0.0 
and we used the Case B rate. 

As noted above, we have defined the cooling rates for 
each species between 5000 and 10^ K. In regions where 
the temperature is below this range, we turned off cooling 
while allowing the chemistry to evolve, and we enforced 
an absolute temperature floor at 100 K. 


4. RESULTS 


Note s. Radia t ive Rec ombi nation ra t es taken 
from iBadnelll (|2006b|> . fllBadnelll 



(iz juiover & Abel| ( |20QH| ' 


The material is assumed to have solar metallicity, and 
each run is defined by an initial uniform density and the 
strength of turbulent forcing. The ionization state of 
each ion is initially set to be consistent with a 10^ K 
gas in collisional ionization equilibrium, except for the 
lowest velocity dispersion runs where all species were as¬ 
sumed to be neutral with an initial temperature of 10^ 
K. For most cases, the eddy turnover timescale was much 
shorter than the timescale for the chemistry to come 
to equilibrium. To prevent long run times, we imple¬ 
mented a method of accelerating the chemistry such that 
it reached the steady state solution in a much shorter 
number of cycles. Once the turbulence had reached a 
steady state, normally after a few eddy turnover times, 
we carried out a ‘kick cycle’ over which the chemistry 
was evolved for a much longer time than either the local 
hydrodynamic time step or the eddy turnover timescale. 
The result of this procedure was to force each cell to 


4.1. Model Parameters 

Our goal is to study the steady-state ionization struc¬ 
ture of a gas that is being stirred with ID velocity dis¬ 
persions between 6 and 58 km s“^ (corresponding to 3D 
velocity dispersions between 10 and 100 km s“^), over a 
wide range of densities. In particular, we are interested 
in cases in which the heating from the turbulent stirring 
is balanced by atomic cooling. 

The parameter space spanned by our simulations is 
greatly simplified by the dependencies of turbulent de¬ 
cay and cooling on density and length scale. For a given 
steady state, the average turbulent energy dissipation 
rate per unit volume, or consequently, the average heat¬ 
ing rate per unit volume is 

T = erg s“^ cm“^, (8) 

where a is the velocity dispersion, p is the average den¬ 
sity, and L is the driving scale of the turbulence. Con¬ 
versely, the average cooling rate per unit volume is 

A oc p‘^X{T)/{pmH)‘^ erg s“^ cm“^, (9) 

where p is the average particle mass and A(T) is the 
average cooling rate. Equating these two terms gives 

A(T) (X ap{Lp), (10) 
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TABLE 2 

Summary of the models. Those models with the appended _C denote runs made with thermal conduction while those 

APPENDED WITH _HIGH DENOTE HIGH RESOLUTION MODELS WITH TWICE THE NORMAL RESOLUTION. THOSE DENOTED WITH AN ASTERISK 

DENOTE MODELS THAT RESULT IN THERMAL RUNAWAY. 


Name 

P 

Column 

(TlD 

Tmw 

Tvw 

Mmw 

Mvw 

s 

(Js 

Sskew 

Skurt 

C’sjcxp 

N1E14_S6* 

7e-31 

l.le+14 

5.8 

5.52 

5.51 

0.29 

0.30 

0.004 

0.05 

-0.60 

1.01 

0.15 

N1E15_S6 

7e-30 

l.le+15 

5.8 

2.55 

2.54 

0.29 

0.30 

0.005 

0.05 

-0.50 

1.74 

0.15 

N1E16_S6 [Low] 

7e-29 

l.le+16 

5.8 

1.49 

1.48 

0.51 

0.51 

0.005 

0.13 

-0.84 

1.90 

0.25 

N1E17_S6 

7e-28 

l.le+17 

5.8 

1.26 

1.24 

0.61 

0.62 

0.001 

0.16 

-0.79 

1.48 

0.30 

N1E16_S12* 

7e-29 

l.le+16 

11.5 

26.85 

26.85 

0.21 

0.21 

0.002 

0.02 

-0.54 

1.95 

0.10 

N3E16_S12 

2e-28 

3.0e+16 

11.5 

7.01 

6.99 

0.44 

0.45 

0.006 

0.08 

-0.60 

1.43 

0.22 

N1E17.S12 

7e-28 

l.le+17 

11.5 

1.30 

1.21 

1.06 

1.15 

-0.061 

0.46 

-0.57 

0.62 

0.53 

N3E17_S12 

2e-27 

3.0e+17 

11.5 

1.14 

1.02 

1.26 

1.39 

-0.119 

0.61 

-0.83 

2.16 

0.63 

N7E16_S20* 

5e-28 

7.6e+16 

20.2 

34.05 

33.99 

0.36 

0.37 

0.006 

0.06 

-0.60 

1.95 

0.18 

N1E17_S20 [Med] 

le-27 

1.5e+17 

20.2 

5.54 

5.31 

0.91 

0.96 

-0.021 

0.31 

-0.68 

1.16 

0.46 

N7E17_S20 

5e-27 

7.6e+17 

20.2 

1.01 

0.82 

2.39 

2.95 

-0.381 

0.98 

-0.19 

-0.27 

1.07 

N1E18_S20 

le-26 

1.5e+18 

20.2 

0.87 

0.69 

2.78 

3.56 

-0.497 

1.11 

-0.21 

-0.11 

1.20 

N4E17_S35* 

3e-27 

4.6e+17 

34.6 

176.12 

175.95 

0.31 

0.31 

0.004 

0.05 

-0.49 

1.34 

0.16 

N1E18_S35 

le-26 

1.5e+18 

34.6 

1.12 

0.89 

4.13 

5.25 

-0.798 

1.42 

-0.33 

0.11 

1.44 

N4E18_S35 

3e-26 

4.6e+18 

34.6 

0.90 

0.73 

4.63 

6.49 

-0.877 

1.52 

-0.37 

0.16 

1.56 

N1E19_S35 

le-25 

1.5e+19 

34.6 

0.83 

0.64 

5.32 

7.65 

-0.844 

1.43 

-0.17 

-0.10 

1.66 

N3E18_S58* 

2e-26 

3.0e+18 

57.7 

67.55 

64.53 

0.91 

0.96 

-0.005 

0.32 

-0.54 

0.73 

0.46 

N1E19_S58 [High] 

7e-26 

l.le+19 

57.7 

1.06 

1.09 

7.98 

9.83 

-1.196 

1.77 

-0.32 

0.07 

1.80 

N3E19_S58 

2e-25 

3.0e+19 

57.7 

0.86 

0.71 

8.79 

13.55 

-1.375 

1.93 

-0.32 

-0.16 

1.96 

N1E20_S58 

7e-25 

l.le+20 

57.7 

0.77 

0.64 

10.64 

14.51 

-1.441 

2.01 

-0.46 

0.33 

2.00 

NlE16_S6_High 

7e-29 

l.le+16 

5.8 

1.49 

1.48 

0.45 

0.45 

0.009 

0.10 

-0.44 

1.12 

0.22 

N1E16_S6_C 

7e-29 

l.le+16 

5.8 

1.50 

1.49 

0.49 

0.49 

0.007 

0.12 

-1.11 

3.35 

0.24 

NlE17_S20_High 

le-27 

1.5e+17 

20.2 

5.66 

5.43 

0.89 

0.94 

-0.021 

0.32 

-0.68 

0.90 

0.45 

N1E17_S20_C 

le-27 

1.5e+17 

20.2 

5.61 

5.35 

0.88 

0.93 

-0.028 

0.34 

-0.64 

0.84 

0.44 

NlE19_S58_High 

7e-26 

l.le+19 

57.7 

0.95 

0.81 

7.68 

10.64 

-1.173 

1.74 

-0.23 

-0.22 

1.84 

N1E19_S58_C 

7e-26 

l.le+19 

57.7 

1.07 

1.14 

7.36 

10.41 

-1.303 

1.94 

-0.57 

0.66 

1.83 


Notes, p is the mean density in units of gm cm Column is the column density in units of cm gid is the 1-D 
velocity dispersion in units of km/s. Tmw and Tvw are the mass-weighted and volume-weighted temperatures in units 
of 10^ K. Mmw and Mvw are the mass-weighted and volume-weighted turbulent Mach numbers, s is the volumed 
averaged value of s = In p/p. cr^, Sskew and Skurt are the rms, skewness and kurtosis excess of the volume-weighted 
probably density function of s. crs,exp is the expected variance from Eqn. 


or in other words that the overall thermal distribution 
of the gas is only dependent on a and the column den¬ 
sity. For the velocity dispersions we are interested in, 
for which cooling balances heating, this column density 
ranges between 10^^ and 10^^ cm“^. 

Furthermore, because collisional ionization and recom¬ 
bination rates are also proportional to density squared, 
the ratio of the eddy turnover time defined here as 
teddy = ^box/2cr, to the dcnsity-temperaturc averaged 
recombination time, tree, is also purely a function of a 
and column density, with tree ^ teddy foi* many species. 
This means that many of the species of interest experi¬ 
ence a wide range of conditions within a recombination 
time, which, as we shall see below in more detail, can 
lead to abundance ratios that can not be predicted from 
local collisional equilibrium. 

A summary of the simulation runs is given in Table 
with the name of each run referring to its column density 
and ID velocity dispersion. In the analysis that follows, 
we have chosen three representative models to describe 
in detail: N1E16_S6, N1E17_S20, and N1E19_S58, which 
span a range of densities and velocity dispersions. We 
will refer to these as Low, Medium, and High respec¬ 
tively, as they represent cases with progressively higher 
velocities and turbulent Mach numbers. Figure shows 


the hydrodynamic and chemical evolution of N1E17_S20 
in units of the eddy turnover timescale. The effect of 
the chemical kick is apparent at teddy ~ 5 as the large 
jump in the fractional abundances of each species, e.g., 
Fc = FeJ After this adjustment, the species 

quickly find a global steady state that is distinct from 
instantaneous collisional ionization equilibrium. To en¬ 
sure that the chemical kick did not introduce any bias in 
our model, we ran model N1E17_S20 without the kick. 
The statistical and atomic properties compared very well 
between these two models, although the run without the 
kick took over four times longer to reach a final steady- 
state. 

Figure shows the density, temperature, and the frac¬ 
tion of carbon in the commonly observed state on 
slices extracted from each of our representative models 
after they have reached a global steady state. Here we see 
that while density and temperature variations are small 
at low turbulent Mach numbers, the more highly tur¬ 
bulent simulations display a wide range of temperatures 
and densities, which are only weakly correlated with each 
other. Similarly, the distribution of the fraction of 
which has a recombination time that is comparable to 
^eddy, displays features that sometimes follow the tem¬ 
perature distribution, but sometimes show substantial 
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Fig. 2. — Hydrodynamical and chemical evolution of N1E17_S20. Top Left: The (red) dash-dotted line shows the maximum temperature, 
the (blue) dashed line shows the minimum temperature, the (black) dotted line shows the volume weighted temperature, and the (black) 
solid shows the density weighted temperature. Top Center: The (blue) solid line shows the minimum density and the (red) dotted line 
shows the maximum density. Top Right: The (blue) solid line shows the total energy, the (red) dotted line shows the internal energy, and 
the (dashed) green line shows the kinetic energy. Middle Row: The fractional abundances of the chemical species where each panel shows 
a different element as noted by the panel legends. Bottom Left: Same as middle row for Oxygen. Bottom Center: Same as the middle row. 
For convenience we have combined Na and Mg. Bottom Right: The (blue) solid, left axis line shows the Mach number of the simulation 
and the (red) dashed, right axis line is the ID velocity dispersion. 


variations, indicating that is often very different 

than would be estimated from local collisional ionization 
equilibrium. Below we study each of these effects in turn. 

4.2. Probability Density Functions 

The top panel of Figure shows the two-dimensional 
volume-weighted probability density functions, which 
quantify the fraction of the volume in the simulations 
located at various temperatures and densities. Consis¬ 
tent with the expectations from Figure the gas in the 
Low simulation has a very small spread of temperatures 
and densities. In the Medium simulation, on the other 
hand, the temperatures and densities span roughly an 
order of magnitude. In this case, for which the (mass 
weighted) turbulent Mach number is Mmw = 0.91, T is 
roughly proportional to as would occur exactly in a 
constant pressure medium, but occurs here with a con¬ 
siderable spread, consistent with the presence of tran¬ 
sonic motions. Finally, in the High case, with a Mach 
number of Mmw = 8.0, the densities and temperatures 


span nearly six orders of magnitude and four orders of 
magnitude, respectively. Furthermore, these quantities 
are largely uncorrelated with each other, as would be ex¬ 
pected in a medium filled with strong, supersonic shocks. 

While our simulations are the first to include full 
rate equations and cooling for a large number of atomic 
species, several studies of isothermal, supersonic turbu¬ 
lence have show n that the gas approxiniates a lognor¬ 
mal distribution (|Vazquez-Semadeni||1994| [Padoan et al. 


|Klessen||2U00l lOstriker et al.ll^^OOll |Li et al.||‘iOO;Il 

Kritsuk et al.| 20071 |Federrath et al.| 2008 

Lemaster &; 

Stone 20081 

ISchmidt et al. 2009| |Federrat 

1 et al.| 2010 

Glover et a 

l.| 120101 IPadoan & Nordlundl 

20111 ICollins 

et al. 20TTf 

Price et al. 2011 Molina et al. 2012p, de- 


fined as 


Ps ds = 




=exp 




(s - SqY 
2(7? 


( 11 ) 


where s is the logarithmic density, s = ln{p/p). 

The mean logarithmic density in this case is related to 
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Fig. 3. — Density (top row), temperature (middle row) and fraction (bottom row) on a fixed slice taken from our three representative 
models, Low (left column), Med (center), and High (right) 


the standard deviation as sq = —cr^/2. The density vari¬ 
ation at a particular location is produced by the succes¬ 
sive passage of shocks with mach numbers independent 
of the local density, which gives a physical explanation 
for the lognormal density distribution. For an isother¬ 
mal distribution, the variance of the logarithmic density 
corresponds to 


= ln(l + 


( 12 ) 


where 6 is a constant t hat depends on the forc ing that 
drives the turbulence. Federrath et al. (2008) showed 
that 6 = 1 for purely compressive, V x F = 0, forcing, 
while b = 1/3 for purely solenoidal, V • F = 0, forcing. 

The bottom, left panel of Figure H shows the proba¬ 
bility density functions of the logarithmic density, which 
take on a gaussian profile for runs Low, Med, and High. 
As expected, the width of the profile increases as the 
steady state Mach number increases, with Low having 
the smallest width and High having the largest. In Ta¬ 
ble we show several important statistical quantities in 
terms of the logarithmic density for each run. In particu¬ 
lar, we calculate the first four moments of the logarithmic 
density distribution: the mean, variance, skewness, and 
kurtosis. 

The mean and variance have the usual definitions. 


while the skewness, {{s — sq)^) measures the sym¬ 

metry of the distribution. A positive skewness denotes a 
distribution that has a long tail toward higher densities 
while a negative skewness represents a long tail toward 
lower densities. The kurtosis, ((s — sq)^) /cr^ quantifies 
the peakedness of the distribution and measures the im¬ 
portance of the tails versus the peak. A larger kurtosis 
represents a flatter distribution while a lower value de¬ 
notes a narrower distribution, with a gaussian distribu¬ 
tion having a kurtosis equal to 3. In Table we give the 
kurtosis excess where this factor of 3 is subtracted. 

Figureshows the variance, skewness, and kurtosis ex¬ 
cess as a mnction of Mach number for each model. The 
solid line in the top panel shows the expected crg-M re¬ 
lation from equation \T^ with b = 0.5. We find that for 
supersonic flows the models match very well with equa¬ 
tion!^ but find a significant mismatch for subsonic flows. 
A least-squares fit of equation to our points suggests 
a good fit to the M > 1 points with b = 0. 53. How¬ 
ever, this do es not match the expectation from |Federrath| 


et al. 


(2008) who suggest a value of 5 = 1/3 for purely 
solenoidal driving, for isothermal turbulent models. This 
difference may simply be due to the non-isothermal na¬ 
ture of our models. The local variation in temperature 
allows for a larger spread in density than in isothermal 
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Fig. 4. — Top: Two-dimensional probability density function for simulations Low (Left), Med (Center), and High (Right). Temperature 
is given along the y-axis in units of K and the number density is given along the x-axis. Here we have normalized by the average number 
density, n, which corresponds to -4.43, -3.02, and -1.94 for Low, Med, and High respectively. All contours are labeled by their values relative 
to the PDF bin with the most mass. Bottom Left: One-dimensional probability density functions for run Low, Med, and High. The ar-axis 
gives the logarithmic density. Low is shown by the (red) solid line, Med is shown by the (blue) dashed line, and High is shown by the 
(green) dotted line. The density variance increases as the stirring increases. Also note that the profiles are not symmetric with longer tails 
toward lower densities. Bottom Right: One-dimensional temperature PDFs for Low, Med, and High. The a^-axis gives the temperature. 
The lines have the same meanings as in the other panel. 


models. 

The right panel of Figure]^ show the 1-D temperature 
PDFs for runs Low, Med, and High. As expected, the 
width of the PDF grows with the strength of the stirring, 
meaning that isothermality is not a good approximation 
for our high mach number runs, which may help explain 
the discrepancy in b. Note also the large temperature 
spread causes a few cells in run High to reach the tem¬ 
perature floor, due to the extreme PdV work done by 
the stirring, but this is only a minor effect. 

The top row of Figure shows the two dimensional 
PDFs of density and temperature for runs Low, Med, 
and High. For the Low case, the slight stirring has only 
produced very minor density contrasts, leaving a mostly 
uniform medium. However, there is a slight tail towards 
lower density and lower temperature. This leads to a 
highly concentrated phase diagram where essentially all 
the mass has the same temperature. This can also be 
seen in the left column of Figure where we show slices 
of density, temperature, and the Sundance of 

A similar trend is seen for run Med in the center panel 
of Figureand the center column of FigureAlthough, 
as mentioned below, the average mach number of the flow 


is higher than Low; turbulence has not produced any 
strong density contrasts. In fact, the density is almost 
always within a factor of 3 of the mean. As in the case 
for Low, the density and temperature slices are largely 
uniform. 

Finally, the High case is quite different than either Low 
or Med. Large density and temperature contrasts are 
seen, due to the much higher final mach number (M = 
8.0). This produces a very wide density-temperature 
PDF as shown in the right panel of Figure 0 with den¬ 
sity and temperature spanning nearly six orders of mag¬ 
nitude and four orders of magnitude respectively. These 
contrasts are plainly seen in the right column of Figure [3 

Finally, the bottom row of Figure shows slices of 
Fc 3+. In general the density and temperature distribu¬ 
tions are very similar. For example, the density and tem¬ 
perature slices for High show that the complex density 
structure is largely mirrored in the temperature. How¬ 
ever, this is not the case for In fact, there is only 

a weak relation between these quantities. As we discuss 
below, this highlights the need to follow the nonequilib- 
rium atomic chemistry exactly when estimating atomic 
ionization states. 
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Fig. 5.— Statistical measures of the s = ln(p/p) PDF for each 
model. The x-axis is the density weighted Mach number. Top: 
The points show the variance in s for eac h m odel, and the solid 
line shows the expected relation from eqn |12| with b = 0.5. Note 
that for supersonic flows the solid line matches very well with the 
points, but for subsonic flows the line does not match the data well. 
Middle: The skewness for each model. Note that for all models, 
the skewness is negative, which indicates longer tails toward lower 
densities. Bottom: The kurtosis. At subsonic Mach numbers most 
models show a more peaked distribution while at supersonic Mach 
numbers the distribution is either nearly gaussian or slightly flatter 
than gaussian. 


4.3. Ionization States 

The large spread in both density and temperature in 
our simulations, particularly in the supersonic cases, nat¬ 
urally translates into a large spread in chemical species. 
In fact, the instantaneous chemical makeup in a given cell 
is a strong function of not only the current temperature 
and density but its thermodynamic history. 

To get an idea of how the global steady state abun¬ 
dances depend on our full treatment of turbulence, cool¬ 
ing, and chemical reactions, we calculated approximate 
abundances using two approaches. For our first method, 
labeled Fiviean, we assumed the full medium was in colli- 
sional ionization equilibrium at the mass weighted aver¬ 
age temperature, and we calculated the ionization states 
for each species using Cloudy. This method quantifies 
the errors involved with assuming constant temperature 
and density conditions in a region with significant tur¬ 
bulent motions. In our second estimate, labeled FrurbEq, 
we used the final temperature PDF to calculate a mass- 
weighted average abundance for each species, again using 
Cloudy to calculate collisional ionization abundances for 
the full range of conditions. This method quantifies the 
errors involved with failing to track the nonequillibrium 
abundances accurately by including reaction rates within 
the simulations. 

The left panel of Figure shows the comparison be¬ 
tween each of these estimates to the steady state abun¬ 
dances for Low, Med, and High. The name of each 
species is listed along the x-axis and the fractional species 
abundances are given along the ^-axis. These values are 
also listed in Table 3. As shown in above, for the Low 
simulation, with a mean temperature of 1.5 x lO^K and 
a Mach number of 0.5, the spread in the temperature- 


density PDF is small. Hence, the differences between 
the global species abundances and those from our two 
estimates are mostly minor, with most of the true val¬ 
ues matching the simpler calculations within a factor of 
2 for most elements with significant abundances. Two 
notable exceptions to this agreement are for which 
the true abundance is ~ 30 times greater than the pMean 
and FrurbEq, values, and He+, which has a long recom¬ 
bination time, and a true abundance ^ 20 times greater 
than expected from more simple estimates. 

In the Med case, for which the mass-weighted mean 
temperature is Tmw = 5.5 x 10^ and M = 0.92, the 
difference between the true values and simple estimates 
becomes more pronounced. Here the true abundances of 
species such as He, C, C+, and Mg, differ by a 

factor of ~ 2 from our simple estimates, while He^+ and 
are more abundant by ~ 5 than expected from either 
approximate approach. Note also that there is no clear 
trend for the more complicated FrurbEq estimate to be 
a better predictor than the simple constant temperature 
approach. 

Finally, in the High case, with Tmw = IT x lO^K 
and M = 8.0, the abundances of many species differ 
by several order of magnitude from simple expectations. 
For example, the abundance of H+ is over two orders of 
magnitude higher than it is in the constant temperature 
models, but close to the TxurbEq estimate. The abun¬ 
dance of He, on the other hand, is very near the value 
expected from TMean and roughly twice the TxurbEq esti¬ 
mate, while and He^+ is many orders of magnitude higher 
than expected by the TMean estimate, but still over 200 
times less abundant than expected in the TxurbEq esti¬ 
mate. Similarly, and 

which are almost completely absent in the TMean, are 
all present at significant levels, although less so than ex¬ 
pected by TxurbEq- On the other hand, Mg+ is present 
at a level much less than expected by the TMean esti¬ 
mate, but more than predicted by TxurbEq- These re¬ 
sults not only show that the presence of turbulence can 
introduce significant abundance differences when M ^ 1, 
and change the abundances completely when M > 1, 
but that these differences can only be predicted by full 
nonequilibrium calculations such as the ones carried out 
here. 

To get an idea of the cause of the large discrepancy 
between carbon and oxygen in the High case, we compare 
the ionization and recombination timescales of these ions. 
This is shown in the right panel of Figure We find that 
while the recombination time is short compared to the 
eddy turnover, the ionization timescale for these ions is 
long. Physically, this suggests that even though a parcel 
of gas is heated due to the crossing shocks, the ions do 
not have sufficient time to fully ionize. This explains the 
lower abundances found for TxurbEq compared to TMean 
for C^+through and through for case High. 

4.4. Thermal Conduetion and Resolution 

In our simulations described above, energy transport is 
purely by advection. Yet even in a supersonic medium, 
the large difference in the thermal velocities of the elec¬ 
trons and the ions can lead to significant energy trans¬ 
port through conduction. In the case in which many col¬ 
lisions take place over the scale length for temperature 
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Fig. 6. — Left: Steady state species fractions for each of the representative runs. Each species is given along the a^-axis, and the y axis 
gives loglO of the fraction of the total mass in a given element that is found in a particular ionization state. The (red) triangles show Fmean, 
the expected species fraction assuming collisional ionization equilibrium at the mass weighted mean temperature. The (blue) diamonds 
show FxurbEqj the fraction assuming collisional ionization equilibrium using the temperature PDF. The (black) stars show, -Epurb? the final 
global steady state fractions. Right: Comparison between the ionization and recombination timescales for each ion. The (red) squares and 
(blue) pentagons are the ratios of the recombination and ionization timescales with the eddy turnover time, respectively. 


variations, the thermal conductivity is given by 
1/2 


K = i^[^] = 5.2xlO“®Ty^n-i 


cm^ s (13) 


where z/ is the plas ma viscosity and T is the temp erature 
in units of 10^ K (Braginskii||l958 Spitzer||l962 ). Com¬ 
paring this to the velocity dispersion and length scale of 
the turbulence gives G\^LjK oc ai^nL. Like the ratios 
of timescales discussed in section 4.1, this is dependent 
only on the column density and turbulent velocity. In 
fact, the Reynolds number of the medium also only de¬ 
pends on these two quantities, although in practice the 
numerical viscosity in our simulations is much greater 
than the true physical value. 

In the case in which the scale length for tempera¬ 
ture variations is smaller than the collisional mean free 
path of the electrons, the conduction becomes saturated. 
The mean free path is = 1.3 x lO^^cm^^Ty 
where rii^c is the ion density, the saturated flux is 

Fsat(^e\e'ne^/{kT)^/iTie^ where (aele is the electron conduc¬ 
tivity flux-limiter coefficient, Ue is the ele ctron number 
density, a nd rup is the mass of the electron (Cowie & Mc- 
Kee|1977 ). From these scalings we can see that L/Xi and 
/nk'l {a) are purely functions of column density and 
the temperature structure of the medium. This means 
that neither unsaturated or saturated conduction intro¬ 


duces a new parameter that must be spanned by our 
simulations. 

To test the impact of electron conduction, we re¬ 
ran our representative simulation using of the implicit 
SpitzerHighZ conductivity module within FLASH, a 
Larsen flux limiter, and an electron conductivity flux- 
limiter coefficient of aeie = 0.2. The upper panel of Fig¬ 
ure compares the logarithmic density PDFs from our 
nominal models, with these conductions runs, from these 
runs, which were roughly twice as expensive in terms of 
computer time. For run Low the PDF profile is essen¬ 
tially unchanged from the nominal model, in particular 
with regards to the mean and variance of s. Table 
presents the statistical measures between our nominal 
runs and those with electron conduction. There is a slight 
increase in the skewness and kurtosis which produces a 
slightly longer distribution to lower densities and slightly 
more peaked distribution. Run Med shows no substan¬ 
tial change in the kurtosis, but conduction does produce 
a slightly less peaked distribution. Run High shows a 
similar trend as the other models, with a slight increase 
in both the skewness and kurtosis, but now with a no¬ 
ticeable difference in the extreme low-density tail of the 
log density PDF. 

The top right panel of Figure shows the effect of elec¬ 
tron conduction on the logarithmic temperature PDF. 
Like the density PDFs, the temperature PDFs are largely 
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TABLE 3 

Species fractions in our representative runs. 



Low 

FMean 

Low 

TpurbEq 

Low 

Ururb 

Medium 

FMean 

Medium 

UrurbEq 

Medium 

Frurb 

High 

FMean 

High 

FfurbEq 

High 

Frurb 

H 

-0.23 

-0.18 

-0.12 

-3.47 

-3.86 

-3.32 

-0.37 

-0.00 

-0.42 

H+ 

-0.39 

-0.46 

-0.63 

0.00 

0.00 

-0.00 

-0.24 

-2.35 

-0.21 

He 

-0.00 

0.00 

0.00 

-2.54 

-2.36 

-1.57 

-0.03 

0.00 

-0.23 

He+ 

-3.26 

-4.44 

-4.49 

-0.07 

-O.OI 

-0.03 

-1.14 

-8.03 

-0.64 

He2+ 

-17.28 

-30.00 

-30.00 

-0.86 

-I.6I 

-1.40 

-3.07 

-30.00 

-0.73 

c 

-0.67 

-0.51 

-0.35 

-3.10 

-3.24 

-2.78 

-0.37 

-0.02 

-0.48 

C+ 

-O.II 

-0.16 

-0.26 

-0.73 

-0.66 

-0.47 

-0.25 

-1.32 

-0.45 

C2+ 

-3.67 

-5.15 

-5.19 

-O.IO 

-O.II 

-0.18 

-1.75 

-9.79 

-0.74 

C3+ 

-14.19 

-30.00 

-18.41 

-1.98 

-2.55 

-2.33 

-3.31 

-30.00 

-1.51 

C4+ 

-18.56 

-30.00 

-30.00 

-3.47 

-5.94 

-4.12 

-4.15 

-30.00 

-0.99 

0 

-0.28 

-0.22 

-0.14 

-3.07 

-2.91 

-2.61 

-0.37 

-0.00 

-0.42 

0+ 

-0.33 

-0.40 

-0.56 

-0.44 

-0.31 

-0.22 

-0.25 

-2.24 

-0.44 

0^+ 

-7.62 

-8.81 

-8.49 

-0.21 

-0.29 

-0.41 

-2.02 

-15.66 

-0.84 

03+ 

-17.97 

-30.00 

-30.00 

-1.88 

-2.68 

-2.38 

-3.20 

-30.00 

-1.10 

04+ 

-18.92 

-30.00 

-30.00 

-5.92 

-30.00 

-6.05 

-4.73 

-30.00 

-1.63 

05+ 

-18.72 

-30.00 

-30.00 

-11.67 

-30.00 

-11.69 

-6.21 

-30.00 

-2.33 

0®+ 

-18.85 

-30.00 

-30.00 

-13.14 

-30.00 

-30.00 

-7.53 

-30.00 

-1.98 

0^+ 

-19.68 

-30.00 

-30.00 

-20.89 

-30.00 

-30.00 

-15.17 

-30.00 

-10.37 

Na 

-4.26 

-4.27 

-4.07 

-5.86 

-6.29 

-6.11 

-2.10 

-3.37 

-1.63 

Na+ 

-0.00 

0.00 

0.00 

-0.00 

0.00 

-0.00 

-0.00 

-0.00 

-0.01 

Mg 

-1.40 

-1.29 

-1.12 

-6.37 

-6.63 

-5.09 

-0.65 

-0.58 

-0.74 

Mg+ 

-0.18 

-0.08 

-0.07 

-3.41 

-3.49 

-2.94 

-0.31 

-0.13 

-0.61 

Mg2+ 

-0.53 

-0.94 

-I.II 

-0.00 

-0.00 

-0.00 

-0.55 

-3.24 

-0.24 


identical for Low and Med. Only for the High simulation 
is there a small decrease in the PDF at the highest tem¬ 
peratures at which conduction is most efficient, which 
is accompanied by an increase in the low-density end 
of the density pdf. Thus it appears that in the highest 
Mach number cases, conduction is able to operate quickly 
enough to move a noticeable amount of energy out of the 
hottest cells, moving some of this energy into low den¬ 
sity regions which are able to expand a bit more due to 
the increase in pressure. However, these differences are 
small and suggest that electron conduction does not play 
a major role. 

In order to simulate a large number of turbulent con¬ 
ditions for many eddy turnover times in the presence 
of chemistry at reasonable CPU cost, our simulations 
have adopted a relatively low resolution of 128^. In this 
case, two-point statistics such as the velocity and density 
power-spectra and structu re functions are not expected 
to be well reproduced (e.g. , ^^ kovichlloM Kritsuk et al 


2007 Pan fc Scannapie^|20I0| , and we do not attempt 


to analyze them here. On the 3ther hand, the probabil¬ 
ity distribution function is a single-point quantity that is 
much easier to simulate at moderate resolution. 

As a check of convergence, the bottom left row of Fig¬ 
ure compares the PDF obtained in our representative 
simulations with a second set of simulations with the 
same set of physical parameters, but a resolution of 256^. 
Again, we find that the resulting PDFs are nearly iden¬ 
tical to the nominal runs and only in the High case is 
there a slight decrease in the high temperature tail and 
a slight dip in the low density tail. However, in terms of 
the mean and variance, the effect is minor, which sug¬ 
gests that our nominal model with N = 128 is sufficient 


to produce the key one-point statistical quantities and 
obtain estimates of the species abundances. 

Finally, in Figure we show the difference in the 
species abundances between our nominal runs and the 
comparison runs including electron conduction and an 
increase in resolution. The comparisons, which are also 
presented in Table show that neither electron conduc¬ 
tion nor an increase in resolution has a substantial effect 
on our abundance results. 


4.5. Parameter Dependeneies 

Having explored the evolution of three representative 
cases in detail, established convergence of abundances at 
a resolution at 128^, and shown that conduction has only 
a minor effect, we turn our attention to the full suite of 
models. Table gives the overall parameters of each of 
these runs, which span a large range of columns, tem¬ 
peratures, and Mach numbers. In most cases, we found 
that heating from stirring and cooling from recombina¬ 
tion are quickly balanced, as seen in the Med simulation 
presented in Fig.[^ However, if the stirring is sufficiently 
strong that the average temperature of the model exceeds 
lO^K these models will undergo thermal runaway. This 
can be explained by the fact that most eleme nts have 
peaks in their cooling functions at ^ IQ^K (e .g.jGnat fc 
Ferland||20I2l jOppenheimer & Schaye| 20I3|) . Therefore 
once the stirring passes this barrier, cooling can no longer 
balance turbulent heating and a meaningful steady state 
can no longer be achieved. 

A summary of our results is given in Table which 
shows the final abundances for all models that were able 
reach a steady state. Inspection of these abundances 
reveals many interesting trends affecting commonly ob- 
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Fig. 7.— Effect of resolution and thermal conduction. Top Left: Effect of electron conduction on the logarithmic density PDF. The 
solid lines are the fiducial models while the dashed lines show the model with conduction. Top Right: Effect of electron conduction on the 
temperature PDF. Bottom Left: Comparison between our fiducial models and models with twice the base resolution on the logarithmic 
density PDF. Bottom Right: Comparison between our fiducial models and models with twice the base resolution on the temperature PDF. 
The solid lines are the fiducial models while the dashed lines show the higher resolution models. 


served species. For example, without any contribution 
from photoionization, He+ is often found at substantial 
levels, even when the average temperature is lO^K, and 
in many of these cases the fraction can exceed 50% 
(as seen when cfid = 35 km/s and Mmw = 2.8, or when 
(JiD — 58 km/s and A^mw = 8.0). 

On the other hand, mostly appears at substantial 
levels when the mean temperature is significantly higher 
than lO^K, such as in the aiD = 11.5 km/s, Tmw = 
7.0 X lO^K, and cjid = 20 km/s, Tmw = 5.5 x lO^K 
runs. Note however, that it is also seen in the high-mach 
number cjuj = 58, km/s run with Tmw = 1.1 x lO^K, 
even though in CIE it is only expected when T > 5 x lO^K 
as shown in Fig. [l] Similarly, which has a slightly 

lower ionization potential, is most abundant in runs with 
elevated temperatures or high mach numbers, although 
at even higher fractions than seen in 

Interestingly, although its ionization potential is only 
5.1 eV, substantial levels of neutral sodium are found in 
most of our simulations, particularly those with higher 
Mach numbers. Inspections plots of the physical distri¬ 
bution of this ion show that it is mostly found in under- 
dense regions, which drop to relatively low temperatures 


through adiabatic cooling in expansions. Finally, Mg+ 
which recombines at 7.6 eV and is ionized at 15 eV is ex¬ 
tremely abundant in all our simulations, with the excep¬ 
tion of the highest temperature aiD = 11.5 km/s Tmw = 
7.0 X lO^K, and aiD = 20 km/s Tmw = 5.5 x lO^K cases. 
In fact, in the absence of photoionization, the existence 
of significant, cospatial Mg+, and is a clear in¬ 
dicator of a broad temperature-density PDF, caused by 
supersonic turbulence. 

At highest Mach numbers, when the dispersions in 
temperature and density are large, regions in the sim¬ 
ulation are able to reach conditions in which the fraction 
of molecular hydrogen may not be negligible. Because 
parcels of gas move between regions with significantly 
different conditions on the timescale in which molecules 
are formed, calculating the H 2 fraction exactly would 
require solving a compete set of rate equations for H 2 
formation within our simulations themselves. While this 
is beyond the scope of this work, we can nevertheless 
make a rough estimate of the molecular fractions calcu¬ 
lating cell by cell abundances of H“, HJ, and H 2 in the 
approximate case in which the local conditions last long 
enough for these additional species to come into equilib¬ 
rium with the species tracked exactly by our simulations. 
In this case, in each cell we can compute 
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TABLE 4 

Effect of resolution and electron conduction of the final elemental abundances. Each row gives a difference species 

WHILE EACH COLUMN GIVES A DIFFERENT SIMULATION. _HiGH DENOTES MODELS WITH TWICE THE BASE RESOLUTION AND THOSE WITH _C 

DENOTE MODELS WITH ELECTRON CONDUCTION. 


Col. (cm ^) 
aiD (km s“^) 

1.1E16 

5.8 

l.lE16_High 

5.8 

1.1E16_C 

5.8 

1.5E17 

20.2 

1.5E17_High 

20.2 

1.5E17_C 

20.2 

1.1E19 

57.7 

l.lE19_High 

57.7 

1.1E19_C 

57.7 

H 

-0.23 

-0.22 

-0.15 

-3.47 

-3.50 

-3.47 

-0.37 

-0.27 

-0.29 

H+ 

-0.39 

-0.40 

-0.52 

0.00 

0.00 

0.00 

-0.24 

-0.34 

-0.31 

He 

-0.00 

-0.00 

-0.00 

-2.54 

-2.59 

-2.56 

-0.03 

-0.02 

-0.05 

He+ 

-3.26 

-3.34 

-3.71 

-0.07 

-0.07 

-0.07 

-1.14 

-1.30 

-0.95 

He2+ 

-17.28 

-17.39 

-17.92 

-0.86 

-0.81 

-0.84 

-3.07 

-3.14 

-2.48 

c 

-0.67 

-0.69 

-0.60 

-3.10 

-3.14 

-3.11 

-0.37 

-0.34 

-0.36 

C+ 

-0.11 

-0.10 

-0.13 

-0.73 

-0.75 

-0.74 

-0.25 

-0.27 

-0.26 

C2+ 

-3.67 

-3.60 

-3.80 

-0.10 

-0.09 

-0.09 

-1.75 

-2.26 

-1.92 

C3+ 

-14.19 

-14.15 

-18.30 

-1.98 

-1.91 

-1.95 

-3.31 

-3.77 

-3.20 

C4+ 

-18.56 

-18.41 

-18.13 

-3.47 

-3.18 

-3.63 

-4.15 

-4.26 

-3.59 

0 

-0.28 

-0.26 

-0.17 

-3.07 

-3.10 

-3.07 

-0.37 

-0.27 

-0.30 

0+ 

-0.33 

-0.35 

-0.50 

-0.44 

-0.46 

-0.44 

-0.25 

-0.34 

-0.31 

02+ 

-7.62 

-8.05 

-13.47 

-0.21 

-0.20 

-0.21 

-2.02 

-2.38 

-1.94 

03+ 

-17.97 

-17.82 

-17.52 

-1.88 

-1.85 

-1.93 

-3.20 

-3.62 

-3.02 

04+ 

-18.92 

-18.73 

-18.44 

-5.92 

-5.85 

-6.01 

-4.73 

-4.81 

-4.54 

05+ 

-18.72 

-18.57 

-18.29 

-11.67 

-11.25 

-11.84 

-6.21 

-5.96 

-6.09 

0®+ 

-18.85 

-18.66 

-18.49 

-13.14 

-10.67 

-16.48 

-7.53 

-7.05 

-6.88 

0"+ 

-19.68 

-19.17 

-19.82 

-20.89 

-19.84 

-20.17 

-15.17 

-17.10 

-17.03 

Na 

-4.26 

-4.29 

-4.29 

-5.86 

-5.85 

-5.85 

-2.10 

-2.12 

-2.22 

Na+ 

-0.00 

-0.00 

-0.00 

-0.00 

-0.00 

-0.00 

-0.00 

-0.00 

-0.00 

Mg 

-1.40 

-1.44 

-1.36 

-6.37 

-6.39 

-6.36 

-0.65 

-0.60 

-0.65 

Mg+ 

-0.18 

-0.17 

-0.12 

-3.41 

-3.44 

-3.41 

-0.31 

-0.24 

-0.24 

Mg2+ 

-0.53 

-0.54 

-0.72 

-0.00 

-0.00 

-0.00 

-0.55 

-0.76 

-0.69 


kiFuF^- +/c23i^H+^e- 

(/c2 + ki^)F}i + {k^ + ^16)^H+ + ki^F^- + (/C21 + ^22)^h+ ^ ^28^He+ + ^29^He 
_ ksFY[FY[+ + k7FY[^FY[+ + + fe5-^H2-^He+ 

2 k^Fji + keF^- + {k2i + k22)Fii- 

^7^H+ + {kg + k2g)FQ- + /CqFh + kiiF}iQ + (/C24 + ^25)^He+ 


(14) 

(15) 

(16) 


where the temperature dependent reaction rates, /ci, 
/c 2 , etc. are taken from the compila tion presented in Ap¬ 
pendix A of Glover & Abel (2008). Furthermore, to es¬ 
timate the formation ol II 2 from neutral-neutral sur face 
reactions on dust we follow Glover & Jappsen (2007), eq. 
(36) which assumes that every collision between a metal 
atom and a grain results in a reaction. In that case, the 
reaction rate per unit volume for an atomic species i is 


given by kR = 3.0 x 10-^^tV2/(1.0 + 4 x + 

2 X 10“^T + 8.0 X 10“^T^). The resulting fractions are 
shown in Table denoted by asterisks because the are 
only approximate measurements. While the majority of 
hydrogen in our simulations remains atomic in all cases, 
these estimates indicate that in the highest Mach number 
runs, ~ 1% of the gas could be in the form of H 2 . 

In addition to Table 0 we make available online the 
data files from each rmjj Each file presents a variety 


^ http://zofia.sese.asu.edu/~evan/turbspecies/ 


of information, including the true abund ances and the 
simple abundance estimates presented in §4.3[ We also 
give Doppler parameters for each species, i, defined as. 


(17) 


where aipD is the ID velocity dispersion, k^ is the Boltz¬ 
mann’s constant, and Ai is the ion atomic mass. While 
the masses of hydrogen and helium are small enough that 
the temperature term makes a substantial contribution, 
the Doppler parameters of the heavier elements are very 
close to cr^pD in general, providing a good measurement 
of the local velocity dispersion. 


5. CONCLUSIONS 

Turbulence pervades the Universe, often moving at su¬ 
personic speeds due to the high efficiency of radiative 
cooling. These random motions provide overall support 
against gravity, but also concentrate a portion of the 
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TABLE 5 

Final steady-state abundances for our suite of models. Abundance values are given as \ ogiQ { Fi / FT ). Here, H“, H+, H2 are 

COMPUTED ACCORDING TO EQS. AND ARE THUS APPROXIMATE. 


Col . (cm 
aiD (km s “^) 
Tdw (lO^K) 
Mdw 

1 . 1 E 15 

5.8 

2.55 

0.29 

1 . 1 E 16 

5.8 

1.49 

0.51 

1 . 1 E 17 

5.8 

1.26 

0.61 

3 . 0 E 16 

11.5 

7.01 

0.44 

1 . 1 E 17 

11.5 

1.30 

1.06 

3 . 0 E 17 

11.5 

1.14 

1.26 

1 . 5 E 17 

20.2 

5.54 

0.91 

7 . 6 E 17 

20.2 

1.01 

2.39 

1 . 5 E 18 

20.2 

0.87 

2.78 

1 . 5 E 18 

34.6 

1.12 

4.13 

4 . 6 E 18 

34.6 

0.90 

4.63 

1 . 5 E 19 

34.6 

0.83 

5.32 

1 . 1 E 19 

57.7 

1.06 

7.98 

3 . 0 E 19 

57.7 

0.86 

8.79 

1 . 1 E 20 

57.7 

0.77 

10.6 

H 

- 1.03 

- 0.23 

- 0.03 

- 3.96 

- 0.19 

- 0.10 

- 3.47 

- 0.17 

- 0.09 

- 0.42 

- 0.19 

- 0.09 

- 0.37 

- 0.15 

- 0.07 

H+ 

- 0.04 

- 0.39 

- 1.13 

0.00 

- 0.45 

- 0.68 

0.00 

- 0.48 

- 0.73 

- 0.21 

- 0.44 

- 0.73 

- 0.24 

- 0.53 

- 0.81 

He 

- 0.00 

- 0.00 

- 0.00 

- 2.87 

- 0.17 

- 0.01 

- 2.54 

- 0.00 

- 0.00 

- 0.01 

- 0.00 

- 0.00 

- 0.03 

- 0.00 

- 0.00 

He+ 

- 2.20 

- 3.26 

- 4.58 

- 0.02 

- 0.49 

- 1.84 

- 0.07 

- 2.96 

- 3.68 

- 1.59 

- 2.72 

- 3.30 

- 1.14 

- 2.61 

- 2.94 

He2+ 

- 15.48 

- 17.28 

- 19.53 

- 1.27 

- 4.52 

- 7.37 

- 0.86 

- 9.96 

- 10.90 

- 5.22 

- 7.76 

- 9.17 

- 3.07 

- 5.66 

- 5.78 

C 

- 1.51 

- 0.67 

- 0.21 

- 3.71 

- 0.67 

- 0.46 

- 3.10 

- 0.50 

- 0.28 

- 0.63 

- 0.34 

- 0.19 

- 0.37 

- 0.19 

- 0.12 

c+ 

- 0.02 

- 0.11 

- 0.41 

- 1.05 

- 0.11 

- 0.19 

- 0.73 

- 0.17 

- 0.33 

- 0.12 

- 0.27 

- 0.44 

- 0.25 

- 0.44 

- 0.61 

C2 + 

- 2.15 

- 3.67 

- 6.06 

- 0.05 

- 3.95 

- 4.54 

- 0.10 

- 4.02 

- 4.50 

- 2.46 

- 3.36 

- 4.03 

- 1.75 

- 3.30 

- 3.75 

C3+ 

- 10.10 

- 14.19 

- 20.35 

- 1.53 

- 12.45 

- 13.21 

- 1.98 

- 11.00 

- 10.54 

- 5.69 

- 7.42 

- 9.28 

- 3.31 

- 6.01 

- 5.65 

C4+ 

- 18.50 

- 18.56 

- 20.01 

- 3.96 

- 16.07 

- 21.04 

- 3.47 

- 14.78 

- 16.02 

- 7.24 

- 10.76 

- 13.39 

- 4.15 

- 6.89 

- 7.15 

0 

- 1.25 

- 0.28 

- 0.05 

- 3.49 

- 0.50 

- 0.25 

- 3.07 

- 0.26 

- 0.10 

- 0.53 

- 0.20 

- 0.09 

- 0.37 

- 0.15 

- 0.08 

0 + 

- 0.03 

- 0.33 

- 0.99 

- 0.64 

- 0.17 

- 0.36 

- 0.44 

- 0.34 

- 0.70 

- 0.15 

- 0.43 

- 0.71 

- 0.25 

- 0.53 

- 0.79 

02+ 

- 4.60 

- 7.62 

- 13.81 

- 0.13 

- 4.39 

- 7.06 

- 0.21 

- 5.74 

- 6.40 

- 3.17 

- 4.58 

- 5.47 

- 2.02 

- 4.20 

- 4.44 

03+ 

- 18.07 

- 17.97 

- 17.73 

- 1.54 

- 14.15 

- 15.07 

- 1.88 

- 11.64 

- 11.91 

- 5.62 

- 8.00 

- 10.23 

- 3.20 

- 5.82 

- 5.70 

04+ 

- 19.46 

- 18.92 

- 18.46 

- 5.51 

- 24.67 

- 25.59 

- 5.92 

- 20.29 

- 19.12 

- 8.90 

- 12.32 

- 16.24 

- 4.73 

- 7.68 

- 7.61 

©5 + 

- 19.54 

- 18.72 

- 18.22 

- 12.31 

- 26.39 

- 26.34 

- 11.67 

- 23.90 

- 26.03 

- 11.26 

- 16.89 

- 21.61 

- 6.21 

- 9.13 

- 10.90 

06 + 

- 19.51 

- 18.85 

- 18.54 

- 18.62 

- 26.33 

- 26.33 

- 13.14 

- 26.34 

- 26.34 

- 14.61 

- 22.27 

- 26.33 

- 7.53 

- 10.96 

- 14.93 

o"+ 

- 21.43 

- 19.68 

- 19.64 

- 20.34 

- 26.54 

- 26.54 

- 20.89 

- 26.54 

- 26.55 

- 26.55 

- 26.56 

- 26.57 

- 15.17 

- 19.93 

- 23.88 

Na 

- 5.31 

- 4.26 

- 3.87 

- 5.76 

- 3.88 

- 3.51 

- 5.86 

- 2.91 

- 2.45 

- 2.75 

- 2.19 

- 1.98 

- 2.10 

- 1.86 

- 1.57 

Na + 

- 0.00 

- 0.00 

- 0.00 

- 0.00 

- 0.00 

- 0.00 

- 0.00 

- 0.00 

- 0.00 

- 0.00 

- 0.00 

- 0.00 

- 0.00 

- 0.01 

- 0.01 

Mg 

- 2.97 

- 1.40 

- 0.98 

- 7.23 

- 1.29 

- 1.04 

- 6.37 

- 0.88 

- 0.67 

- 0.98 

- 0.61 

- 0.48 

- 0.65 

- 0.44 

- 0.35 

Mg+ 

- 0.83 

- 0.18 

- 0.06 

- 3.84 

- 0.16 

- 0.09 

- 3.41 

- 0.13 

- 0.12 

- 0.25 

- 0.18 

- 0.20 

- 0.31 

- 0.25 

- 0.29 

Mg2+ 

- 0.07 

- 0.53 

- 1.57 

- 0.00 

- 0.60 

- 1.03 

- 0.00 

- 0.93 

- 1.42 

- 0.48 

- 1.00 

- 1.38 

- 0.55 

- 1.13 

- 1.46 

H“* 

- 8.39 

- 7.42 

- 7.39 

- 11.79 

- 7.36 

- 7.32 

- 11.13 

- 7.38 

- 7.36 

- 7.61 

- 7.42 

- 7.37 

- 7.58 

- 7.40 

- 7.40 

hU 

- 7.08 

- 6.90 

- 7.51 

- 9.38 

- 6.99 

- 7.14 

- 9.07 

- 7.08 

- 7.20 

- 7.10 

- 7.11 

- 7.22 

- 7.18 

- 7.16 

- 7.30 

H 2 * 

- 7.85 

- 6.24 

- 6.15 

- 14.08 

- 5.84 

- 3.90 

- 12.84 

- 3.16 

- 2.14 

- 2.86 

- 2.17 

- 1.84 

- 2.21 

- 1.86 

- 1.60 


material to very high densities, giving rise to a multi¬ 
phase distribution with unique thermodynamic proper¬ 
ties. This can have a dramatic effect on the chemical 
makeup of the medium. Specifically, if the recombina¬ 
tion time for a given species is long compared to the 
eddy turnover time, it cannot reach an equilibrium state 
before it is further acted upon by the turbulence. This 
creates a situation in which the final ionization state is 
not only a function of the temperature and density, but 
also a function of the rate at which parcels of gas move 
through these conditions. 

To study how these abundances are altered by tur¬ 
bulence, we have implemented a nonequilibrium atomic 
chemistry package within FLASH. This package tracks 
the evolution of six elements and 24 separate ionization 
states. In additi on, we have used the method of |Gnat| 


& Ferland (2012) to derive ion-by-ion cooling curves for 
each of the 10 ns under consideration that allows us to fol¬ 
low the thermodynamic evolution of the gas. The result 
is a very fast and efficient package, which we are able to 
run for many eddy turnover times for many cases. 

Using this tool, we have performed a suite of direct 
numerical simulations of solenoidally-driven turbulence 
over a range from ID velocity dispersion of criD= 6 to 
58 km s“^ and the product of the mean density and 
turbulent length scale from 10^^ to 10^^ cm“^ for so¬ 
lar metallicity gas, concentrating on three representative 
models. As found by isothermal models of driven tur¬ 
bulence, the gas approximates a lognormal distribution, 
whose logarithmic density variance in the supersonic case 
is well approximated by = ln(l + 6^M^). On the other 
hand, this expression overestimates the variance at sub¬ 
sonic Mach numbers, and the 6 = 0.53 value that best fits 
our data is significantly different from the 5 = 1/3 value 
measured in solenoidally-driven, isothermal turbulence. 


We compare the final steady state abundances in our 
simulations to those obtained assuming the gas is in colli- 
sional ionization equilibrium, using both the mean tem¬ 
perature and the full temperature PDF. We find that 
at low mach numbers the estimates agree to within a 
factor of ^ 2 for most species, save for He^+ and 
which show large deviations due to their long recombi¬ 
nation times. At intermediate mach numbers, several 
species such as He, C, C+, and Mg, differ by 

a factor of ~ 2 from our simple estimates, while He^+ 
and are ~ 5 more abundant than our simples esti¬ 
mates suggest. Finally, for very high mach numbers, the 
abundances can vary by many orders of magnitude from 
simple estimates. Neither increasing the resolution by a 
factor of two or including of thermal electron conduction 
has a significant effect on these abundances. 

These results underscore the fact that transsonic and 
supersonic turbulence can drastically alter the abun¬ 
dances and that only nonequilibrium calculations can 
predict these changes accurately. Thus we make make 
all of the derived properties from our models available 
online. In particular, we present the logarithmic den¬ 
sity statistics and other hydrodynamic quantities, such 
as Mach number and average temperature. We also give 
the final abundances for each species, the abundance val¬ 
ues from the two simple estimates, and species by species 
Doppler parameters. 

In future work we, plan on increasing our network to 
include ions of nitrogen (N-N^+), silicon (Si-Si^+), and 
iron (Fe-Fe^+) and their associated cooling as well as the 
effects of photoionization. Using this increased network, 
we will run similar models as those presented here, again 
compiling important statistical properties and the final 
ionization structure of the gas. This will result in a large 
set of tables, useful for a variety of theoretical and ob- 
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Fig. 8.— Effect of electron conduction and resolution. Each 
species is given along the x-axis. The (red) triangles show the 
steady state species fractions with electron conduction. The (blue) 
diamonds show the steady state species fractions with twice the 
base resolution. The (black) stars show the final global steady 
state. 


servational applications. 
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